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Reaction-diffusion systems may lead to the formation of steady state heteroge- 
neous spatial patterns, known as Turing patterns. Their mathematical formulation 
is important for the study of pattern formation in general and play central roles in 
many fields of biology, such as ecology and morphogenesis. In the present study 
we focus on the role of Turing patterns in describing the abundance distribution of 
predator and prey species distributed in patches in a scale free network structure. 
We extend the original model proposed by Nakao and Mikhailov by considering food 
chains with several interacting pairs of preys and predators. We identify patterns of 
species distribution displaying high degrees of apparent competition driven by Tur- 
ing instabilities. Our results provide further indication that differences in abundance 
distribution among patches may be, at least in part, due to self organized Turing 
patterns, and not necessarily to intrinsic environmental heterogeneity. 
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I. INTRODUCTION 

Reaction-diffusion systems, in which two or more species interact locally and diffuse 
through the medium, have long been focus of studies in many different fields, such as Physics, 
Chemistry and Biology. Part of the interest in these systems is related to their potential to 
form self-organized spatio-temporal patterns, like traveling and spiral waves pQ or stationary 
patterns, called Turing patterns [2] . 

Working on the problem of morphogenesis [3] Turing derived general analytical conditions 
for the formation of stationary patterns in reaction-diffusion systems under a mechanism 
today called diffusion driven instability (or Turing instability). Despite the specific nature of 
the original problem, the work led to a large number of applications in chemistry and biology, 
both theoretical [1HS] and, more recently, empirical [TUHT2] . A key theoretical contribution 
was provided by Mimura and Murray [5J, who applied Turing's idea to understand patchiness 
in continuously distributed predator-prey populations. 

Recently, Nakao and Mikhailov JT3] proposed a discrete version of the prey-predator 
model of Mimura and Murray [6J in which the species are organized in patches, instead of 
being continuously distributed in space. The patches are represented by nodes of a complex 
network such that predators and preys interact locally in each patch and diffusion occurs 
through connected nodes. The Turing patterns obtained in [13] present significant differences 
when compared to the ones obtained in the analogous system which considers space as a 
continuous medium. 

In the present work, we extend of the model of Nakao and Mikhailov [T3] by considering 
food chains with more than two species. We study the dynamics of several pairs of preys and 
predators that interact by consuming common preys. We show that the Turing patterns of 
population density displayed by the system present nontrivial correlations in the abundance 
distributions. In particular, we observe the emergence of strong competition between preys 
of adjacent species in the food chain, despite the fact that no direct competition between 
them are included in the equations. These correlations are strictly related to diffusion and 
correspond to a new mechanism of apparent competition, driven by Turing instabilities 
instead of local interactions. We characterize these patterns using numerical simulations 
and mean field approximations. We also discuss the relevance of these results to patterns of 
species distribution in real trophic systems. 
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II. DYNAMICAL MODEL 



In order to consider more complex reaction diffusion systems we extend the model in- 
troduced by Nakao and Mikhailov [13] to food chains composed by several species of preys 
and predators. We assume that each prey species has a primary predator associated to it, 
forming a pair. The pairs in the food chain are hierarchically coupled by secondary preda- 
tion relations. Thus, the prey in the first pair is consumed by its main predator and also 
by the predator in the second pair, though with the lower intensity 7. Similarly, the prey of 
the second pair is consumed primarily by its associated predator and also by the predator 
of the third pair, and so on. Only the last species of prey in this ordered chain is consumed 
exclusively by its main predator as, illustrated by the diagram in figure [T] 

The environment where these interactions take place consists of a network of patches. 
Species-species interactions, as described by the food chain, occur locally in each patch and 
the coupling between patches is exclusively due to diffusion, which is possible if the patches 
are connected in the network. 

preys predators 




FIG. 1. (Color online) Hierarchical food chain with three pairs of preys and predators. Each 
predator is linked to the previous prey (secondary predation) with strength 7. 

The equations describing this dynamical system are given by: 

j 

jvf{t) = g(^,v?) + + <re (1) 

j 
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where uf\t) and vf\t) represent the populations of preys and predators at time t respec- 
tively. The label I = 1, 2, 3, • • • index the prey-predator pair. 

The functions / and g describe the local interaction between preys and predators of each 
pair (also called reaction functions) . The terms proportional to 7 represent the secondary 
predation relations between adjacent pairs and the parameter accounts for the ratio be- 
tween predator gain and prey loss in the secondary interaction. 

The parameters e and a are, respectively, the prey mobility and the ratio between predator 
and prey mobilities. The matrix L stands for the Laplacian matrix and accounts for the 
diffusion of populations across connected sites. For undirected networks L is symmetric 
with Ly = Aij — kiSij, where A is the Adjacency Matrix and ki the degree of node i. The 
adjacency matrix defines the topology of the network and is given by A^ = 1 if nodes i 
and j are connected and A^ = if they are not. The degree ki = J2j A43 is the number of 
connections of node i. 

The term ^2jLijUj in eq, 1 controls the diffusion of preys u®. It gives the difference 
between the total population of preys in the sites connected to i and ki times the 
population in the site i. If is the same in all sites the sum adds to zero and there is no 
diffusion. A similar term controls the diffusion of predators in the equation for v"\ 

As a simplification, we consider that the intrinsic growth rate of all prey species are the 
same, as is the intrinsic death rate of all predator species. In that manner, the functions / 
and g and the parameters associated to these functions are the same for all pairs. 

The functions / and g are chosen according to the model of Mimura and Murray [B] : 



,'a + bu-u 2 , 
f{u,v) = [ v ) u 



(2) 



g(u,v) — [u — (1 + dv)]v, 
where a, b, c and d are positive parameters that will be fixed to a = 35, b = 16, c = 9 and 
d = 0.4 throughout this paper [6]. 

Both the prey per capita growth rate and the pradator per capita death rate are density 
dependent. The hump effect that can be noted in the prey growth in / represents what in 
Biology is called the Allee Effect [T9TI2T] . describing a positive correlation between population 
density and per capita growth rate in small populations. The linear function related to 
the predator per capita death rate accounts for intraspecific competition in the predator 
population. 
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The possibility of observing Turing patterns for these equations must be evaluated via 
linear analysis. Here we show the analysis for the case of a single prey-predator pair. The 
general case with n pairs is slightly more complicated, but can be done following the same 
lines. 



III. LINEAR STABILITY ANALYSIS 

In this section we briefly review the stability analysis of network organized systems. For 
simplicity we consider only one pair of prey and predator, since the methodology generalizes 
immediately to the case of multiple pairs. 

The equilibrium populations in the absence of diffusion, (u,v), are the positive solution 

of: 

/(M = o (3) 

g(u,v) = 

For diffusion driven instability to take place, the equilibrium must be stable against 
small perturbations in the absence of diffusion (e = 0.0) and go unstable, when diffusion is 
considered. 

Let 

(Ui,Vi) = (u, v) + (dui, 5vi) (4) 

be small perturbations to the fixed point (u, v) at site i. Substituting Q in Q and lineariz- 
ing, we obtain 

N 

dt 



^-8ui = f u 5ui + f v 5vi + e ^ L i:j 5ui 



— 5vi = g u 5ui + g v 5vi + ere ^ Lijdvi 

where the derivatives are evaluated at the equilibrium. 

Since we are dealing with network-organized systems, it is convenient to expand the 
perturbations in the basis formed by the eigenvectors of the Laplacian matrix, {$ Q } [13J. 
Here a = 1, ■ ■ • , N represent different modes, in direct analogy with the Fourier modes that 
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appear in continuous systems where the Laplacian is the usual operator V 2 . We find 

TV 

5ui(t) = c a exp[X a t)<p. 

a=l 
N 

5vi(t) = c a B a exp[\ a t](f), 



N 

(,(<*) 

(6) 



a=l 



Substituting (6) in (5) and using Ylf=i ^v^j = A*$i , we obtain, for each mode a: 

y y y ^ + aeA a J yB a J 
The matrix obtained in ^ is the Jacobian of the system with diffusion. The linear growth 
rates, X a , of each mode are, as expected, the eigenvalues of the Jacobian matrix. Turing 
instability appears when one of the modes becomes unstable. At the threshold, Re(X a ) = 
for some a = a c and Re(X a ) < for all other modes. 

Above this threshold Re(X ac ) > and perturbations grow in time according to exp[X a t], 
eventually forming the stationary Turing pattern. A necessary condition for this is that the 
solutions of ([5]) are confined, otherwise the perturbation would diverge. 

Figure [2] shows the linear growth rates, X a , as a function of the eigenvalues of the Jacobian, 
Aq,, when we consider the functions (|2]), with parameters a = 35.0, b = 16.0, c = 9.0, 
d = 0.4 and e = 0.06, and a network of N = 200 nodes with power law degree distribution 
constructed according to the Barabasi- Albert algorithm [H]. Below the critical value o c = 
15.5 (see appendix [A]) X a < 0.0 for all the modes and the homogeneous state (|3| is stable. 

IV. 2 SPECIES 

We first review the case of two species as a reference to the more complex patterns we 
study in the following sections. We consider a network with iV = 1000 nodes, constructed 
according to the Barabasi-Albert model [H]. The populations of preys, Ui, and predators, Vi, 
defined in each node i, interact locally and diffuse through the network nodes according to 
the equations (jl]), with I = 1 (and uf^ = vf^ 1 = 0). Equations (jl| are numerically integrated 
until a stationary distribution of the species abundance is obtained. 



Figure shows the stationary abundance patterns of preys, figure 3(a), and predators, 



figure 3(b)| as a function of node index i, for e = 0.12 and a = 20.0. The nodes are ordered 



according to decreasing degree fcj. 
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FIG. 2. (Color online) Linear growth rates, X a , as a function of the eigenvalues of the Laplacian, 
A Q , for a Barabasi- Albert network with N = 200 and (k) = 10. In all the cases e = 0.06 and three 
different values of a are shown for comparison. Modes with A Q > 0.0 are observed for a > a c = 15.5. 
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(a) (b) 

FIG. 3. (Color online) Stationary abundance patterns for a single predator-prey pair as a function 
of node index % for e = 0.12 and a = 20.0 for (a) preys and (b) predators. The lines in Uj = 5.0 
and Vi = 10.0 indicate the values of the homogeneous state, which is a fixed point for this set of 
parameters. 

The pattern of prey distribution is formed by two groups of nodes presenting significant 
differentiation in relation to the homogeneous state: a group with high abundance (values 
of Ui well above u) and a group with low abundance (values of «j well below u). The pattern 
of predators follows directly the pattern of the preys: nodes with large abundance of preys 
(lii > u) also have large abundance of predators (i>j > v) and vice versa. 
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V. 4 SPECIES 

The 4 species system is described by Eq.|l| with I — 1,2 (and uf^ 1 = i>® = 0). The 
equations have a homogeneous equilibrium point that depends on the coupling parameter 7, 
as displayed by the table [TJ The populations of preys and predators decrease as 7 increases. 
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u (l) u (2) v (2) 


0.002 
0.01 
0.05 


4.989 9.973 4.993 9.995 
4.945 9.863 4.966 9.977 
4.726 9.314 4.837 9.889 



TABLE I. Homogeneous fixed points for different values of 7 for the four species system. 




FIG. 4. (Color online) Stationary abundance patterns for it 1 ) and vS 2 ' as a function of node index 
i for e = 0.12, a = 20.0, 4> = 0.5 and 7 = 0.002. 

The stationary patterns of preys and as a function of the node index i are shown 
in figure [ZJ These patterns of abundance (and also those of and v^) are not very 
different from each other or from the previous case shown in figure [3] In particular, both 
types of preys and predators present the separation of nodes in high abundance and low 
abundance groups. 

However, this similarity is partly an illusion, having to do with the way the data is plotted. 
Indeed, a new underlying pattern arises when difference between the prey abundances u'p — 
is plotted, as shown in figured for different values of the coupling strength 7. 
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(a) (b) (c) 

FIG. 5. Stationary patterns for the difference — as a function of node index i for (a) 
7 = 0.002, (b) 7 = 0.01 and (c) 7 = 0.05. In all cases e = 0.12, a = 20.0 and = 0.5. 

In all cases it is possible to distinguish three main branches: the upper branch, where 
u w — u ( 2 ) ~ 4^ corresponding to nodes where is abundant but is not; the lower 
branch, where u^ 1 ' — u^ 2 ' ~ —4 where the abundances are reversed; and the middle branch, 
where vP^ — vP^ ~ and and have similar abundances. This configuration of 
branches can be derived via a mean field approximation [22j [23] , as discussed in appendix 
[B] and displayed in Fig. [6j 
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FIG. 6. (Color online) Stationary pattern for the difference w 1 ' —u^ as obtained from simulations 
(black squares) and from the mean field approximation (blue) for e = 0.12, a = 20.0, <fi = 0.5 and 
7 = 0.05. The cyan lines show unstable branches. 

As 7 increases the middle branch gets less populated and the nodes are dominated mostly 
by a single species of prey and predator. This corresponds to a strong effect of apparent 
competition driven by Turing instabilities. The more important is the secondary predation 
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(which is kept weaker than the direct predation in the each pair), the stronger is the effect. 



VI. 6 SPECIES 



To investigate if the negative correlation between preys of coupled pairs also occur in 
larger trophic chains we consider a system with 6 species, again given by equation ([!]) with 
I = 1, 2, 3 (and uf^ = = 0). The stationary patterns of preys distributions are shown in 
figure [7} 
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FIG. 7. (Color online) Stationary patterns of preys distributions as a function of node index i for 
the case of three pairs for e = 0.12, a = 20.0, (ft = 0.5 and 7 = 0.05. 

Once again, for each prey species, the nodes cluster into groups of high and low abun- 
dances. The analysis of the correlations between different prey species, however, is now more 
involved. We first define the quantity: 

a y = sng{u\> -« w = < , (8) 

[ -1, if uf < tZ« 

where af^ indicates if the Z-th prey population at node i has high (af' = +1) or low 
(crf' ) = — 1) abundance with respect to the homogeneous value. 

Second, we separate the nodes in two groups: those with = +1 and those with 
o~( 2 ) = — 1. Since nodes with large ki are not sensitive to the coupling, we restrict this 
analysis to nodes with i > 250, for which the differentiation is more evident. Finally we 
focus on the value of the sum + for these nodes. The three possible values of this 
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FIG. 8. (Color online) Histogram of nodes with different values of + for = — 1 (black 
bars) and = +1 (red bars) for: (a) 7 = 0.0, (b) 7 = 0.002, (c) 7 = 0.01 and (d) 7 = 0.05 

sum indicate the following situations: if er^ + = +2, both u^ 1 ' and u^ 3 ' have high 
abundance in the node; if + = —2, both w 1 ' and have low abundance and if 
+ = 0, and vS 3 ' have opposed abundance characteristics. If the hypothesis of 
negative correlation is to be valid, the group of nodes with = +1 must have most of its 
node with + = —2 and the group with = —1 must have most of its nodes with 
o-M + or( 3 ) = +2. The results are shown in figure [8] in the form of histograms. 

In Fig. 8(a) 7 = 0.0 and the three preys distributions are uncorrelated. Figures 8(b)| and 



8(c) display cases with increasing values of 7. As the coupling strength increases, the number 
of nodes with + = +2 increases in the group with = — 1 and similarly with the 
number of nodes with + = —2 in the group where = +1. This separation is 



evident in figure 8(d), where 7 = 0.05, where it is clear that most of the nodes where 
has large abundance display low abundances of both and and vice versa, showing 
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the persistence of the negative correlation between preys of coupled pairs. 

VII. DISCUSSION 

We have studied the formation of Turing patterns in an extended prey-predator system, 
considering trophic chains composed of 1, 2 and 3 prey-predator pairs, coupled by cross 
predation and dispersing through the connected nodes of a complex network. We detected 
the emergence of negative correlations between the populations of preys of coupled pairs in 
each node, even though there are no direct competition between preys in the equations. This 
effect, known in Biology as apparent competition [211 EH], is triggered here by the Turing 
instabilities, and not by the local interactions. 

The description of fragmented landscapes as complex networks is relatively recent in 
ecology [15J. Although large landscape networks have been studied [26], most of the empirical 
work has dealt with a relatively small number of patches [27] and it is not obvious that the 
patterns observed here for networks with iV = 1000 nodes persist in smaller sets. We have 
checked that for iV as low as 100 the same pattern of apparent competition can be clearly 
identified, but not so much for N = 50, which seems to be a limiting size for the present set 
of parameters. 

Another important concern in the application of our results to realist ecological problems 
is the topology of the network. All numerical simulations presented in the previous sections 
were performed for networks exhibiting power law decay of the degree distribution, that 
results from the application of the Barabasi-Albert algorithm. Natural landscape networks 
can exhibit significant heterogeneity in the degree distribution [16], but are not necessarily 
scale free. In order to verify the robustness of our results against changes in the network 
topology we have also simulated networks with Poisson degree distribution, associated to 
random networks. We found that the negative correlations between preys still holds for 
iV = 1000 and average degree (k) = 20. 

The occurrence of Turing patterns in real ecological systems is still an open question. 
This is in part due to the difficulties in conducting controlled ecological field experiments 
to distinguish between patterns related to space heterogeneity or to intrinsic mechanisms of 
the interaction. However, there is growing evidence of species distribution patterns formed 
by the Turing mechanism [T7J [18] . 
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Our results point to the possibility that, at least in part, species abundance patterns might 
be related to Turing instabilities and not to environmental heterogeneity. Moreover, strong 
effects of apparent competition might emerge spontaneously as Turing patterns, resulting 
from diffusion instabilities and not necessarily from local interactions. 
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Appendix A: Critical value for Turing instability 

The eigenvalues of the Jacobian matrix of [7] are given by the roots of the characteristic 
polinomial 

X 2 a - X a (fu + g v + (1 + cr)eA a ) + (f u + eA a )(g v + aeA a ) - f v g u = 
which are given by 

, U + 9v + (l + o-)eA a ± y/4f v g u + (f u -g v + (l- a)eA a ) 2 

A a — ^ - { ' 

For each mode a there are two possible values for X a , but only the one associated to 
the plus sign can become positive, so we only need to consider this eigenvalue. Solving 



d(X a )/d(A a ) = we obtain the critical Laplacian eigenvalue. Substituting this value in Al 



and imposing that Re(X ac ) = in the instability threshold, we obtain the critical value a c : 

fu9v — %fv9u + 2a/ fv9u(fv9u — fu9v) /An , 

o- c = y 2 • (A2) 

J u 
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Appendix B: Mean field approximation 

The mean field approximation consists in averaging the heterogeneous degree distribution 
of the network by adjusting the strength by which each node senses the presence of its 
neighbors. Introducing the local fields 



T i _ s^n A m 

x i — 2^=1 A %3 U j 

A..JQ 
Vi — 2^7=1 A \ 



fBll 



HJ "j 



and substituting in Eq.Q, we obtain 



d 

dt 
d 

dt 



uf(t) = S ^\,f)-WPv^ +C { X \-h^) 

«'"(*) = 9(«?\«!") + 07«f~ 1, <f + ^ (»! " ■ (B2) 



We then consider the approximations £■ ~ fcjX and ?/• ~ kiY , where the global fields X 
and Y are defined as the weighted averages 

y i ) (B3) 



with the weights 



v 

»i = *i/I>- ( B4 ) 



This choice gives hubs have a stronger influence in the calculation of the global fields. 

With this approximation, and introducing the parameter f3(i) = eki, the dynamical 
system may be written as: 



= ^(u (,) ,« (0 ) + (fovP-^v® + a(3 (Y l - u«) , (B5) 

where each dynamical variable interacts only with its associated global field. Since every 
node now possesses the same dynamical equation, we may drop the index i. 
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In order to describe the patterns for the difference of prey populations, in the case with 
2 prey-predator pairs, we define the new variables: 



u± 
v± 



u« ± u (2) 

v m± v w 



(B6) 



The system of equations (B5), written with the new variables (B6), is given by: 



— = F±(u-,u + ,v-,v+) - -(«_ +u + )(v+ -v-)+0 (X ± -u+) 



dv± 

— „^ v „_,^,„_,^ / ^ 4 



G ± (u.,u + ,v.,v + ) ± ^(u. + u+)(v + - v.) + a(3 (Y ± - v+) , (B7) 



where 



and 



F ± = f(uU(u^u + ),v( 1 \v^,v + ))±f(uW(u-,u + ),vW(v-,v + )) 
G± = g(uW(u..,u + ), vW(v_, v+)) ± g(u^(u^u + ), v^(v^v + )) 



(B8) 



(B9) 



X ± = X 1 ± X 2 
Y ± = Y 1 ± Y 2 

If the global fields for each dynamical variable are given, the parameter /3 may be seen 
as a bifurcation parameter. It is possible to note a saddle node bifurcation in the system, 
and the appearance of new stable fixed points, when the value of (3 is increased from = 0. 

We obtain the global fields (B9) by numerically integrating equations ([T]) and using the 



stationary values of the dynamical variables in (B3) and these in (B9). We then construct 



bifurcation diagrams calculating, for each value of (3, the fixed points of the system (B7). 
Since each node has an associated degree k iy and, therefore, an associated (3, it is possible 
to project the bifurcation diagram in the stationary pattern that resulted of the numerical 
integration of 0. The projection of the bifurcation diagram relative to the variable w_ on 
the stationary pattern for the difference — is shown in figure [6] 



